.........................................   Wahadło magnetyczne

Autorem projektu jest Tomasz Urbańczyk - student V roku fizyki na Wydziale Fizyki, Astronomii i Informatyki Stosowanej Uniwersytetu Jagiellońskiego

Opis projektu

Celem projektu jest symulacja ruchu wahadła magnetycznego. Wahadło magnetyczne to metalowa kulka wiszącą na nitce (zwykłe wahadło), która dodatkowo porusza się w polu magnetycznym wytworzonym przez układ kilku magnesów  (kulka waha się ponad płaszczyzną na której są zamocowane magnesy). Z fizycznego punktu widzenia mamy tu do czynienia z ruchem oscylatora harmonicznego (tłumionego) znajdującego się w dodatkowym polu sił.
W moim projekcie założyłem, że ruch wahadła odbywa się na płaszczyźnie oddalonej o odległość d od płaszczyzny z magnesami. Wszystkie wielkości fizyczne podane są w jednostkach bezwymiarowych. W mojej pracy opiarłem się na wzorach fizycznych dostępnych na poniższych stronach internetowych:

,,http://www.inf.ethz.ch/personal/muellren/pendulum/index.html#simulation", 
  ,,http://www.codeproject.com/KB/recipes/MagneticPendulum.aspx".

Część I

Na poniższym wykresie dynamicznym można zobaczyć trajektorie wahadła (dla  z góry ustalonego czasu ruchu) w zależności od prostych parametrów. Położenie początkowe wybieramy poprzez kliknięcie dowolnego punktu na płaszczyźnie wykresu lub przeciąganie szarej kropki. Czerwone punkty symbolizują magnesy.

Clear["Global`*"];
Manipulate[
x1=-Sqrt[3]/2;
y1=-0.5;
x2=Sqrt[3]/2;
y2=-0.5;
x3=0;
y3=1;
{x0,y0}=a;
magn=Graphics [{PointSize[Large],Red,Point[{x1,y1}],Point[{x2,y2}],Point[{x3,y3}]}];
p=NDSolve[{x''[t]+R*x'[t]+Cst*x[t]==(x1-x[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(x2-x[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(x3-x[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,
y''[t]+R*y'[t]+Cst*y[t]==(y1-y[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(y2-y[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(y3-y[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,x'[0]==0,x[0]==x0,y'[0]==0,y[0]==y0},{x,y},{t,100}];


cc=ParametricPlot[Evaluate[{x[t],y[t]}/.p],{t,0,100},PlotRange->pr,AxesOrigin->{0,0}];Show[cc,magn],{{a,{0,0}},Locator},{{R,0.2,"stała tłumienia"},0.1,1,0.1},{{d,0.4,"odleglość wahadła od płaszczyzny"},0.2,2,0.2},{{Cst,1,"siła oscylatora"},0.1,2,0.1},{{pr,1.1,"zakres rysunku"},1.1,4.1,0.5}]

Część II

W części tej tworzony jest dynamiczny wykres umożliwiający zmianę (oprócz parametrów obecnych w części I) położenia magnesów, a także obserwację zmiany położenia wahadła w czasie. Wychylenie początkowe wybiera się klikając w dowolny punkt na wykresie lub przeciągając szarą kropkę.

Clear["Global`*"];
Manipulate[
{x1,y1}=n;
{x2,y2}=m;
{x3,y3}=q;
{x0,y0}=a;
magn=Graphics [{PointSize[Large],Red,Point[n],Blue,Point[m],Green,Point[q]}];
p=NDSolve[{x''[t]+R*x'[t]+Cst*x[t]==(x1-x[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(x2-x[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(x3-x[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,
y''[t]+R*y'[t]+Cst*y[t]==(y1-y[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(y2-y[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(y3-y[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,x'[0]==0,x[0]==x0,y'[0]==0,y[0]==y0},{x,y},{t,100}];


cc=ParametricPlot[Evaluate[{x[t],y[t]}/.p],{t,0,tim},PlotRange->pr,AxesOrigin->{0,0}];Show[cc,magn],{{a,{0,0}},Locator},{{R,0.2,"stała tłumienia"},0.1,1,0.1},{{d,0.3,"odleglość wahadła od płaszczyzny"},0.2,2,0.1},{{Cst,1,"siła oscylatora"},0.1,2,0.1},{{pr,1.1,"zakres rysunku"},1.1,4.1,0.5},{{tim,50,"czas"},0.1,50},{{n,{0,0.5},"Magnes czerwony"},{-2,-2},{2,2}},{{m,{-Sqrt[3]/4,-1/4},"Magnes niebieski"},{-2,-2},{2,2}},{{q,{Sqrt[3]/4,-1/4},"Magnes zielony"},{-2,-2},{2,2}},Button["Reset położenia magnesów",{n={0,0.5},m={-Sqrt[3]/4,-1/4},q={Sqrt[3]/4,-1/4}}],Button["Reset parametrow wahadła",{R=0.2,Cst=1,d=0.4,pr=1.1,tim=50}]]

Część III

Ta część pracy umożliwia stworzenie wykresu zależności położenia końcowego wahadła (tzn. zbadania koło którego magnesu zatrzyma się wahadło) w zależności od położenia startowego. Stworzony rysunek jest eksportowany do pliku "rysunek.bmp".  Ilość punktów, dla których przeprowadza się obliczenia, ustawia się zmieniając ręcznie parametr ,,rozdzielczosc" (rysunek ma ,,rozdzielczość" ma x ma punktów).   Uwaga, obliczenia trwają dość długo.

Clear["Global`*"];

(*Parametry Startowe*)
rozdzielczosc=30;
x1=-Sqrt[3]/2;
y1=-0.5;
x2=Sqrt[3]/2;
y2=-0.5;
x3=0;
y3=1;
d=0.4;
R=0.2;
cc=1;

(*inicjalizacja zmiennych*)
maxli=rozdzielczosc^2;
tabred=Table[0,{i,maxli},{j,2}];
tabblue=Table[0,{i,maxli},{j,2}];
tabgreen=Table[0,{i,maxli},{j,2}];
kr=1;
kg=1;
kb=1;

Timing[{
For[li=1,li<maxli,{
x0=-2+4*Quotient[li,rozdzielczosc]/rozdzielczosc;y0=-2+4*Mod[li,rozdzielczosc]/rozdzielczosc;
p=NDSolve[{x''[t]+R*x'[t]+cc*x[t]==(x1-x[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(x2-x[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(x3-x[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,
y''[t]+R*y'[t]+cc*y[t]==(y1-y[t])/(Sqrt[(x1-x[t])^2+(y1-y[t])^2+d^2])^3+(y2-y[t])/(Sqrt[(x2-x[t])^2+(y2-y[t])^2+d^2])^3+(y3-y[t])/(Sqrt[(x3-x[t])^2+(y3-y[t])^2+d^2])^3,x'[0]==0,x[0]==x0,y'[0]==0,y[0]==y0},{x,y},{t,1000}];
zz=1000;
w=Flatten[Evaluate[{x[zz],y[zz]}/.p]];
If[Last[w]>0,{tabred[[kr,1]]=x0,tabred[[kr,2]]=y0,kr++},{If[First[w]<0,{tabgreen[[kg,1]]=x0,tabgreen[[kg,2]]=y0,kg++},{tabblue[[kb,1]]=x0,tabblue[[kb,2]]=y0,kb++}]}];
li++;}]
}]

graf=Graphics[{Red,PointSize[0.02],Point[tabred],Blue,PointSize[0.02],Point[tabblue],Green,PointSize[0.02],Point[tabgreen],Black,PointSize[0.05],Point[{0,0}],Black,PointSize[0.05],Point[{x2,y2}],Black,PointSize[0.05],Point[{x3,y3}],Black,PointSize[0.05],Point[{x1,y1}]}];
osie=Plot[y=0,{x,-2,2},PlotRange->{-2,2}, AspectRatio->1];
Show[{osie,graf}]
Export["obrazek.bmp",%]